{ Image Macros for intensity calculations on Photometrics CCD images }

{Globals}
var
    {16 to 8 conversion info}
    ymin, ymax,
    
    {pid numbers for images}
    customLUTpid,
    raw16Pid, rawWidth, rawHeight,
    proc16Pid, proc8Pid, procXmin, procXmax,
    temp16Pid,
    seg8aPid, seg8bPid,
    darkData16Pid, unif16Pid, darkUnif16Pid,
    
    {YES/NO parameters}
    skipFlat, {nonzero = skip flat field correction}
    skipMedian, {nonzero = skip median filter step}
    skipSmooth, {nonzero = skip smooth step}

    {smooth kernel parameters}
    {flat field ratio constant}
    
    
: integer;
{initialize/restore globals}
begin
    requiresUser('Pixel16u',2);
    requiresUser('GetPutPixel',1);
    requiresUser('timer',1);
    requiresUser('utilities',1);
    requiresUser('markup',1);

    ymin := getMemo('ymin');
    ymax := getMemo('ymax');

    {pid numbers}

    customLUTpid := getMemo('customLUTpid');
    raw16Pid := getMemo('raw16Pid');
    rawWidth := getMemo('rawWidth');
    rawHeight := getMemo('rawHeight');
    proc16Pid := getMemo('proc16Pid');
    temp16Pid := getMemo('temp16Pid');
    proc8Pid := getMemo('proc8Pid');
    procXmin := getMemo('procXmin');
    procXmax := getMemo('procXmax');
    seg8aPid := getMemo('seg8aPid');
    seg8bPid := getMemo('seg8bPid');
    
    { := getMemo('');}
    SetBackgroundColor(0);
    SetForeGroundColor(255);
end;


{ Procedure: }
{   create uniform image from a series of data images. }
{   import dark image. }
{   import data image. }
{   flat field correction. }
{   reduce noise and smooth. }
{   convert to 8 bit with automatic scaling. }
{   adjust 8 bit scaling. }
{   threshold }
{   use wand tool to select a series of fragments. }
{   after each wand click, invoke macro to define the segment }
{     (fills roi on segment image with segment number)}
{   and a separate macro to assign the segment to a class }
{     (places class number into array indexed by segment number)}
{     (different macro for each class)}
{   use wand tool to select any nearby bright spots (dirt) }
{   invoke macro to define dirt segment }
{   circular dilate segments. }
{   copy the segments to another image, convert to one value, dilate more, }
{     subtract to produce background segment definition }
{   measure each segment in 16 bit data for sum, area, standard deviation, min, max}
{   Use "analyze particles" command on 8 bit segment image to find area, coordinates, seg number}
{   output to spreadsheet format text window: }
{     image name, top left coordinates, area, sum, standard deviation. }
{ also need way to identify values for internal standard}


{status window: next segment number}

procedure forceROIWithin;
var
    left, top, rwidth, rheight, iwidth, iheight: integer;
begin
    GetPicSize(iwidth,iheight);
    GetRoi(left,top,rwidth,rheight);
    if rwidth = 0 then selectAll; {this fixes most cases}
    GetRoi(left,top,rwidth,rheight);
    if (left < 0) 
    or (top < 0) 
    or (left + rwidth > iwidth) 
    or (top + rheight > iheight) then begin
        putmessage('ROI must not extend outside image');
        exit;
        {make ROI doesn't hack it if ROI wasn't rectangular...}
        if left < 0 then begin
            rwidth := rwidth + left;
            left := 0;
        end;
        if top < 0 then begin
            rheight := rheight + top;
            top := 0;
        end;
        if left + rwidth > iwidth then begin
            rwidth := iwidth - left;
        end;
        if top + rheight > iheight then begin
            rheight := iheight - top;
        end;
        makeroi(left,top,rwidth,rheight);
    end;
end;

procedure forceUncalib;
begin
    choosePic(proc8Pid);
    if Calibrated then begin
        selectAll;
        copy;
        dispose;
        SaveState;
        setNewSize(rawWidth, rawHeight);
        makeNewWindow('Processed 8 bit image');
        proc8Pid := pidNumber;
        setMemo('proc8Pid',proc8Pid);
        RestoreState;
        Paste;
        KillRoi;
    end;
end;

{adjust xmin/xmax using mean ± stdev}
procedure enhanceStdev;
var
    mean, sigma, coef: real;
begin
    choosePic(proc16Pid);
    KillROI;
    choosePic(proc8Pid);
    forceUncalib;
    KillROI;
    RestoreROI;
    forceROIWithin;
    KillROI;
    coef := (procXmax - procXmin + 1) / (ymax - ymin + 1);
{might not work if coef < 0???}
    linLUT16uto8(customLUTPid, procXmin, procXmax, ymin, ymax);
    Cnvrt16uto8(proc16Pid, customLUTPid, proc8Pid);
    RestoreRoi; {take mean & stdev over ROI of 8 bit image}
    SaveState;
    SetOptions('Area,Mean,Std. Dev.,User1,User2');
    Measure;
    RestoreState; {does not restore option settings???}
    mean := (rmean[rCount]-ymin) * coef + procXmin + coef / 2;
    sigma := rStdDev[rCount] * coef + coef / 2;
    ruser1[rCount] := coef;
    ruser2[rCount] := mean;
    {serious round off errors happen when sigma < coef }
    {so that the mean is not known well enough, }
    {image comes out white or black}
    if sigma < coef then sigma := coef;
    {SetCounter(rCount - 1);}
    procXmin := mean - 2*sigma; {this needs to be an adjustable parameter}
    procXmax := mean + 4*sigma;
    SetMemo('procXmin',procXmin);
    SetMemo('procXmax',procXmax);
end;

{display 16 bit data into the 8 bit image using specified xmin/xmax}
procedure show16;
begin
    choosePic(proc16Pid);
    KillROI;
    choosePic(proc8Pid);
    KillROI;
    linLUT16uto8(customLUTPid, procXmin, procXmax, ymin, ymax);
    Cnvrt16uto8(proc16Pid, customLUTPid, proc8Pid);
    RestoreRoi;
    SelectPic(proc8Pid);
    ShowMessage(procXmin,' min\',procXmax,' max');
end;

macro 'Set 8 bit display range';
begin
    ymin := getnumber('gray level for smallest pixel value',ymin);
    ymax := getnumber('gray level for largest pixel value',ymax);
    SetMemo('ymin',ymin);
    SetMemo('ymax',ymax);
end;

procedure swapTemp16;
var
    temp: integer;
begin
    temp := temp16Pid;
    proc16Pid := temp16Pid;
    temp16Pid := temp;
    choosePic(proc16Pid);
    SetPicName('Processed 16 bit image');
    choosePic(temp16Pid);
    SetPicName('Temporary 16 bit image');
end;

procedure hide8image;
begin
    selectPic(proc8Pid);
    setforegroundcolor(255);
    setbackgroundcolor(0);
    clear;
    writeln('Press 8 to display image');
end;

{if the scratch windows are wrong size or missing, create them}
procedure makeScratchIfNeed;
var
    width, height: integer;
begin
    saveState;
    if (ymin = 0) and (ymax = 0) then begin
        ymin := 254;
        ymax := 1;
    end;
    if (ymin < 0) or (ymin > 255) then ymin := 254;
    if (ymax < 0) or (ymax > 255) then ymax := 1;
    {linLUT16uto8 doesn't work if ymin > ymax (bug?)}
    {(want ymin = 254 and ymax = 1 for low numbers = black)}
    {also the mean±stdev macro code fails if ymin > ymax}
    if ymin > ymax then begin
        ymin := 1;
        ymax := 254;
    end;
    SetMemo('ymin',ymin);
    SetMemo('ymax',ymax);
    if not pidExists(customlutPid) then begin
        setNewSize(256,256);
        makeNewWindow('custom LUT');
        SelectAll;
        KillRoi;
        customLUTpid := pidNumber;
        SetMemo('customLUTpid',customLUTpid);
    end;
    linLUT16uto8(customLUTpid, 0, 65535, ymin, ymax);
    if not pidExists(raw16Pid) then begin
        putMessage('makeScratch no raw16');
        exit;
    end;
    choosePic(raw16Pid);
    getPicSize(width, height);
    rawWidth := (width div 4) * 2;
    rawHeight := height;
    setMemo('rawWidth',rawWidth);
    setMemo('rawHeight',rawHeight);
    if rawWidth * 2 <> width then begin
        putMessage('makeScratch raw width not multiple of 4');
        exit;
    end;

    if pidExists(proc16Pid) then begin
        choosePic(proc16Pid);
        getPicSize(width, height);
        if (width <> rawWidth * 2) or (height <> rawHeight) then begin
            dispose;
        end;
    end;
    if not pidExists(proc16Pid) then begin
        setNewSize(rawWidth * 2, rawHeight);
        makeNewWindow('Processed 16 bit image');
        SelectAll;
        KillRoi;
        proc16Pid := pidNumber;
        SetMemo('proc16Pid',proc16Pid);
    end;

    if pidExists(temp16Pid) then begin
        choosePic(temp16Pid);
        getPicSize(width, height);
        if (width <> rawWidth * 2) or (height <> rawHeight) then begin
            dispose;
        end;
    end;
    if not pidExists(temp16Pid) then begin
        setNewSize(rawWidth * 2, rawHeight);
        makeNewWindow('Temporary 16 bit image');
        SelectAll;
        KillRoi;
        temp16Pid := pidNumber;
        SetMemo('proc16Pid',proc16Pid);
    end;

    if pidExists(proc8Pid) then begin
        choosePic(proc8Pid);
        getPicSize(width, height);
        if (width <> rawWidth) or (height <> rawHeight) then begin
            dispose;
        end;
    end;
    if not pidExists(proc8Pid) then begin
        setNewSize(rawWidth, rawHeight);
        makeNewWindow('Processed 8 bit image');
        SelectAll;
        KillRoi;
        proc8Pid := pidNumber;
        setMemo('proc8Pid',proc8Pid);
    end;

    if pidExists(seg8aPid) then begin
        choosePic(seg8aPid);
        getPicSize(width, height);
        if (width <> rawWidth) or (height <> rawHeight) then begin
            dispose;
        end;
    end;
    if not pidExists(seg8aPid) then begin
        setNewSize(rawWidth, rawHeight);
        makeNewWindow('Segments A');
        SelectAll;
        KillRoi;
        seg8aPid := pidNumber;
        setMemo('seg8aPid',seg8aPid);
    end;

    if pidExists(seg8bPid) then begin
        choosePic(seg8bPid);
        getPicSize(width, height);
        if (width <> rawWidth) or (height <> rawHeight) then begin
            dispose;
        end;
    end;
    if not pidExists(seg8bPid) then begin
        setNewSize(rawWidth, rawHeight);
        makeNewWindow('Segments B');
        SelectAll;
        KillRoi;
        seg8bPid := pidNumber;
        setMemo('seg8bPid',seg8bPid);
    end;

    restoreState;
end;

procedure copyRawToProc;
begin
    choosePic(raw16Pid);
    selectAll;
    copy;
    killRoi;
    choosePic(proc16Pid);
    selectAll;
    paste;
    killRoi;
end;

{if there is no raw data image, import one}
procedure importIfNeed;
var 
   origPid: integer;
begin
    if not pidExists(raw16Pid) then begin
        SaveState;
        SetImport('8-bits,Custom');
        SetCustom(2634,1034,2124);
        Import('');
        origPid := pidNumber;
        {MakeNewWindow will not make odd width windows.}
        {Therefore, 16 bit images must be even # pixels wide}
        {or width multiple of 4}
        SetNewSize(2632,1032);
        MakeRoi(0, 2, 2632,1032);
        Copy;
        MakeNewWindow(GetPicName);
        raw16Pid := pidNumber;
        SetMemo('raw16Pid',raw16Pid);
        Paste;
        KillROI;
        ChoosePic(origPid);
        Dispose;
        RestoreState;
    end;
    makeScratchIfNeed;
    if origPid <> 0 then begin
        CopyRawToProc;
    end;
end;

macro '[1] use proc as flat field dark image for data';
macro '[2] use proc as flat field dark image for uniform';
macro '[3] use proc as flat field uniform image';
macro 'include median filter step';
macro 'skip median filter step';
macro 'include smoothing step';
macro 'skip smoothing step';
macro 'include flat field step';
macro 'skip flat field step';

macro 'start over from raw image';
begin
    importIfNeed;
    CopyRawToProc;
    hide8Image;
    writeln('raw data');
end;

macro '[z] undo last 16 bit transform';
begin
    swapTemp16;
    hide8Image;
    writeln('undo');
end;

macro '[r]reduce noise';
begin
    choosePic(proc16Pid);
    makeRoi(2,1,(rawWidth-2)*2,rawHeight-2);
    choosePic(temp16Pid);
    makeRoi(2,1,(rawWidth-2)*2,rawHeight-2);
    median16u(proc16Pid,temp16Pid);
    choosePic(proc16Pid);
    killROI;
    choosePic(temp16Pid);
    killROI;
    swapTemp16;
    hide8Image;
    writeln('reduce noise');
end;

macro '[m]min spatial filter';
begin
    choosePic(proc16Pid);
    makeRoi(2,1,(rawWidth-2)*2,rawHeight-2);
    choosePic(temp16Pid);
    makeRoi(2,1,(rawWidth-2)*2,rawHeight-2);
    minspat16u(proc16Pid,temp16Pid);
    choosePic(proc16Pid);
    killROI;
    choosePic(temp16Pid);
    killROI;
    swapTemp16;
    hide8Image;
    writeln('min spatial');
end;

macro '[*]Convert to 8 bit with min max scaling';
begin
    importIfNeed;
    minmax16u(proc16Pid, procXmin, procXmax);
    SetMemo('procXmin',procXmin);
    SetMemo('procXmax',procXmax);
    show16;
end;

macro '[8]Convert to 8 bit with mean ± stdev scaling';
begin
    importIfNeed;
    minmax16u(proc16Pid, procXmin, procXmax);
    SetMemo('procXmin',procXmin);
    SetMemo('procXmax',procXmax);
    show16;
    enhanceStdev;
    show16;
end;

macro '[„]Enhance ROI of 8 bit image';
begin
    importIfNeed;
    enhanceStdev;
    show16;
end;

macro '[9]reduce xmin';
begin
    importIfNeed;
    procXmin := round(procXmin - 0.1*(procXmax - procXmin) - 1);
    if procXmin > procXmax then procXmax := procXmin + 1;
    SetMemo('procXmin',procXmin);
    show16;
end;
macro '[»]increase xmin';
begin
    importIfNeed;
    procXmin := round(procXmin + 0.1*(procXmax - procXmin) + 1);
    if procXmin > procXmax then procXmax := procXmin + 1;
    SetMemo('procXmin',procXmin);
    show16;
end;

macro '[0]reduce xmax';
begin
    importIfNeed;
    procXmax := round(procXmax - 0.1*(procXmax - procXmin) - 1);
    if procXmax < procXmin then procXmin := procXmax - 1;
    SetMemo('procXmin',procXmin);
    show16;
end;
macro '[¼]increase xmax';
begin
    importIfNeed;
    procXmax := round(procXmax + 0.1*(procXmax - procXmin) + 1);
    if procXmin > procXmax then procXmin := procXmax - 1;
    SetMemo('procXmin',procXmin);
    show16;
end;

macro 'remove calibration on 8 bit image';
begin
    forceUncalib;
end;

macro 'Show calibration numbers';
begin
    ShowMessage('Analyze/optionCalibrate... straight line',
        '\(hold option key while selecting Calibrate)',
        '\measured ',ymin,' known ',procXmin,
        '\measured ',ymax,' known ',procXmax);
    setCounter(2);
end;

macro 'rename front image';
begin
    SetPicName(GetString('new image name',GetPicName));
end;

macro 'Hilight marked areas in sequence';
var
    fg, i: integer;
begin
    fg := pidNumber;
    ChoosePic(seg8aPid);
    killRoi;
    measure
    for i := 1 to 254 do begin
        if histogram[i] <> 0 then begin
            setDensitySlice(i,i);
            showMessage('Mark number ',i:0);
            KillDelay(1);
            StartDelay(1,1.0);
            WaitDelay(1);
        end;
    end;
    setDensitySlice(0,0);
    SelectPic(fg);
end;


procedure appendROI;
var 
    fg, lower,upper: integer;
begin
    fg := pidNumber;
    GetThreshold(lower,upper);
    SetDensitySlice(0,0);
    KillRoi;
    RestoreRoi;
    Clear;
    ChoosePic(seg8aPid);
    SetBackgroundColor(0);
    SetForegroundColor(fore);
    RestoreRoi;
    fill;
    SetForegroundColor(fore);
    SelectPic(fg);
    ShowMessage('Foreground color is ',fore);
    SetDensitySlice(lower,upper);
end;


macro 'ROI is next segment';
begin
    fore := fore + 1;
    if fore > 250 then
        fore := 1;
    if fore < 1 then
        fore := 1;
    SetMemo('fore', fore);
    appendROI;
end;

macro 'append ROI to segment';
begin
    appendROI;
end;